(* ::Package:: *)

(* Particularidades:
- Codificaci\[OAcute]n de Gray 
- Kappa funci\[OAcute]n c\[UAcute]bica de epsilon1 (en tantos por mil)
- Todas las deformaciones est\[AAcute]n en tantos por mil
- C\[AAcute]lculo de la soluci\[OAcute]n del sistema con c\[AAcute]lculo previo de klim para calcular emax y comprobar la hip\[OAcute]tesis de primera cedencia
- Elitismo, sin preselecci\[OAcute]n de la generaci\[OAcute]n 0
- Desecha soluciones fuera de rango (epsilon1 < ectm o kappa > klim2)
*)
Clear["Global`*"];
(* Environment variables *) 
messages=OpenWrite["c:\\MODEL\\messages.txt"];inputdata="c:\\MODEL\\DataBase_Kuchma_et_al_nonprestressed_mathematica.xls";
validhyp="c:\\MODEL\\Critical_Outlet_kappa_refined.xls";
initpop="c:\\MODEL\\initial_population.csv";
thetaseeds="c:\\MODEL\\Strut-angle_seeds.xls";
epsilon1seeds=3; (* Number of equispaced seeds for epsilon1 to try *)
penal=10^7;    (* Fitness penalization for specimens with no solution *)
outputxls="c:\\MODEL\\Salida.csv";
outputpng="c:\\MODEL\\Salida.png";
specimens={"A75","C50","C75","S1-4","S2-3"};
(*specimens={"NHW-3b","A50","A75","C50","C75","S1-4","S2-3","S4-6","ET3","P20","T-22","S8 A","H 50/4","H 75/4","H 100/4","RC 30 A1","RC 30 A2","RC 60 B1","MHB 2.5-25","T3","T4","T6","T7","T8","T9","T13","T15","T17","T19","T20","T26","T32","T34","T35","T37","T38"}; *)


(* Functions used by the algorithm *)


Population[n_,size_]:=Module[{i,j,P,R},
Init=Import[initpop];
P=Init[[1]];
While[Length[P]<n,
R={};
For[j=1,j<=size,j=j+1,
R=Append[R,RandomInteger[{0,1}]];
];
P=Append[P,R];
];
P=Take[P,n];
Return[P];
];


FPSelection[P_]:=Module[{fit,sum,i,Q,r,j,sump,elite},
fit={};
For[i=1,i<=Length[P],i=i+1,
(* Print["      Computing fitness for member ",i,": ",P[[i]]]; *)
fit=Append[fit,Fitness[P[[i]]]];
];
sum=Sum[fit[[i]],{i,Length[P]}];
(* Write[messages,"Fitnesses: "<>ToString[fit]]; *)
Write[messages,"   Mean Fitness: "<>ToString[sum/Length[P]]<>" ("<>ToString[sum*100/Length[P]]<>"%)   Max. Fitness: "<>ToString[ Max[fit]]<>" ("<>ToString[Max[fit]*100]<>"%)"];
meanfit=Append[meanfit,sum*100/Length[P]];
maxfit=Append[maxfit,Max[fit]*100];
Q={};
elite=Max[fit];
i=0;
While[i<Length[P] && Q=={}  ,
i=i+1;
If[fit[[i]]==elite,Q=Append[Q,P[[i]]]];
];
For[i=Length[Q]+1,i<=Length[P],i=i+1,
r=RandomReal[{0,sum}];
j=0;
sump=0;
While[sump<r,
j=j+1;
sump=sump+fit[[j]];
];
Q=Append[Q,P[[j]]];
];
Return[Q];
];


Crossover[P_,pco_]:=Module[{Q,i,j,size,Ha,Hb,Ta,Tb},
size=Length[P[[1]]];
Q={P[[1]]};
For[i=3,i<=pco*Length[P],i=i+2,
j=RandomInteger[{1,size-1}];
Ha=Take[P[[i-1]],j];
Hb=Take[P[[i]],j];
Ta=Take[P[[i-1]],j-size];
Tb=Take[P[[i]],j-size];
Q=Append[Q,Join[Ha,Tb]];
Q=Append[Q,Join[Hb,Ta]];
];
Q=Join[Q,Take[P,Length[Q]-Length[P]]];
Return[Q];
];
Mutation[P_,pmut_]:=Module[{Q,size,i,irandom,j},
size=Length[P[[1]]];
Q=P;
For[i=1,i<=pmut*Length[Q],i=i+1,
irandom=RandomInteger[{2,Length[Q]}];
j=RandomInteger[{1,size}];
Q[[irandom,j]]=1-Q[[irandom,j]];
];
Return[Q];
];


Fitness[R_]:=1/(1+F[R]);
F[R_]:=Module[{i,j,T,Es,EqSystem,suma,\[Alpha],\[Alpha]1,\[Alpha]2,\[Alpha]t,z,bw,Asx1,Asx2,Ast,s,fyx1,fyx2,fyt,fc,\[Epsilon]c,fctm,\[Epsilon]ctm,Acx1,Acx2,Act,Ec,V,\[Sigma]stexp,sol,\[Kappa]limx1,\[Kappa]limx2,\[Kappa]limt,\[Kappa]lim2x1,\[Kappa]lim2x2,\[Kappa]lim2t,\[Kappa]lim1,\[Kappa]lim2,epsilonx,epsilont,\[Epsilon]yx1,\[Epsilon]yx2,\[Epsilon]yt,check,hypx1,hypx2,hypt,emaxx1,emaxx2,emaxt,solx1,solx2,solt,iter,y0,x0,x1,x2,a,b,c,d,H,spec,line,hyp,h,diff,specdiff,seedvector,S,lineseed,Sol},
a=Coef[B2Real[G2B[ Take[R,{1,Length[R]/4}]]]];
b=Coef[B2Real[G2B[ Take[R,{Length[R]/4+1,2*Length[R]/4}]]]];
c=Coef[B2Real[G2B[ Take[R,{2*Length[R]/4+1,3*Length[R]/4}]]]];
d=Coef[B2Real[G2B[ Take[R,{3*Length[R]/4+1,Length[R]}]]]];
T=Import[inputdata];
Es=200000;
suma=0;
For[spec=1,spec<=Length[specimens],spec++,
line=0;
For[i=1,i<=Length[T[[1]]] && line==0, i++,
If[T[[1,i,1]]==specimens[[spec]],line=i];
]; 
If[line==0 ,
Write[messages,"Specimen "<>ToString[specimens[[spec]]]<>" is not in DataBase. Check if 'path' for DataBase file is right."];
Exit[];
];
\[Alpha]=T[[1,line,7]];
\[Alpha]1=T[[1,line,4]];
\[Alpha]2=T[[1,line,5]];
\[Alpha]t=T[[1,line,6]];
z=T[[1,line,8]];
bw=T[[1,line,9]];
Asx1=T[[1,line,10]];
Asx2=T[[1,line,11]];
Ast=T[[1,line,12]];
s=T[[1,line,13]];
fyx1=T[[1,line,14]];
fyx2=T[[1,line,15]];
fyt=T[[1,line,16]];
fc=T[[1,line,17]];
Ec=8500Power[fc+8, (3)^-1];
\[Epsilon]c=T[[1,line,18]];
fctm=T[[1,line,19]];
\[Epsilon]ctm=fctm/Ec;
Acx1=T[[1,line,20]];
Acx2=T[[1,line,21]];
Act=T[[1,line,22]];
Es=200000;
V=T[[1,line,23]];
\[Epsilon]yx1= fyx1/Es;
\[Epsilon]yx2=fyx2/Es;
\[Epsilon]yt=fyt/Es;
\[Sigma]stexp=T[[1,line,2]];
(* Limits to 'Degradation Parameter' for each type of reinforcement *)
\[Kappa]limx1=Asx1*fyx1/(\[Alpha]1*Acx1*fctm);
If[\[Alpha]2!=0,\[Kappa]limx2=Asx2*fyx2/(\[Alpha]2*Acx2*fctm)];
\[Kappa]limt=Ast*fyt/(\[Alpha]t*Act*fctm);
\[Kappa]lim1=Min[\[Kappa]limt,If[\[Alpha]2!=0,Min[\[Kappa]limx1,\[Kappa]limx2],\[Kappa]limx1]];
\[Kappa]lim2x1=\[Kappa]limx1*(4500*\[Epsilon]yx1+(1+1500*\[Epsilon]yx1)^(3/2)-1)/(6750*\[Epsilon]yx1);
If[\[Alpha]2!=0,\[Kappa]lim2x2=\[Kappa]limx2*(4500*\[Epsilon]yx2+(1+1500*\[Epsilon]yx2)^(3/2)-1)/(6750*\[Epsilon]yx2)];
\[Kappa]lim2t=\[Kappa]limt*(4500*\[Epsilon]yt+(1+1500*\[Epsilon]yt)^(3/2)-1)/(6750*\[Epsilon]yt);
\[Kappa]lim2=Min[\[Kappa]lim2t,If[\[Alpha]2!=0,Min[\[Kappa]lim2x1,\[Kappa]lim2x2],\[Kappa]lim2x1]];
(* Seeds for theta *)
S=Import[thetaseeds];
lineseed=0;
For[i=1,i<=Length[S[[1]]] && lineseed==0, i++,
If[S[[1,i,1]]==specimens[[spec]],
lineseed=i;
seedvector={};
];
];
If[lineseed==0 ,
Write[messages,"Specimen "<>ToString[specimens[[spec]]]<>" has no theta seeds. Check if 'path' for theta seeds file is right."]; 
seedvector={30 Degree,30 Degree,30 Degree,30 Degree,30 Degree};
];
(* Valid behaviour hypotheses for reinforcement  *)
H=Import[validhyp];
line=0;
For[i=1,i<=Length[H[[1]]] && line==0, i++,
If[H[[1,i,1]]==specimens[[spec]],line=i];
]; 
If[line==0 ,
Write[messages,"Specimen "<>specimens[[spec]]<>" has no valid behaviour hypotheses. Check if 'path' for hypotheses file is right."];
Exit[];
];
hyp={};
For[h=0,h<=4,h++,
If[H[[1,line,2+9*h]]!=""  &&  H[[1,line,2+9*h]]!="Null",
hyp=Append[hyp,Table[StringTake[H[[1,line,2+9*h]],{k}],{k,1,3}]];
If[lineseed!=0,seedvector=Append[seedvector,S[[1,lineseed,2+h]]]];
];
];
(* Print["             Hypothesis ",hyp,"   Semillas ",seedvector]; *)
diff={};
For[i=1,i<=Length[hyp],i++,
(* Print["                         ",hyp[[i]],"   Semilla ",seedvector[[i]]]; *)
j=1;
specdiff={};
While[j<=epsilon1seeds, (*** Search for solutions with several seeds for epsilon1 ***)
(* Print["                       ",hyp[[i]]," - j:",j,"  \[Epsilon]1seed=",\[Epsilon]ctm+(Min[\[Epsilon]yx1,\[Epsilon]yx2]+\[Epsilon]yt-\[Epsilon]ctm)*(j-1)/(epsilon1seeds-1),"  \[Theta]seed=",seedvector[[i]]]; *)
Clear[\[Sigma]1,\[Sigma]2,\[Sigma]sx1,\[Sigma]sx2,\[Sigma]st,\[Epsilon]1,\[Epsilon]2,\[Epsilon]x,\[Epsilon]t,\[Theta],f2max,\[Kappa]];
\[Sigma]1=\[Alpha]*fctm/(1+Sqrt[500*\[Epsilon]1]);
\[Sigma]2=(Tan[\[Theta]]+1/Tan[\[Theta]])*V/(z*bw)-\[Sigma]1;
f2max=Min[fc,fc/(0.8+170*\[Epsilon]1)];
\[Epsilon]2=\[Epsilon]c*(1-Sqrt[1-\[Sigma]2/f2max]);
\[Epsilon]t=(\[Epsilon]2*Tan[\[Theta]]^2+\[Epsilon]1)/(Tan[\[Theta]]^2+1);
\[Epsilon]x=\[Epsilon]1+\[Epsilon]2-\[Epsilon]t;
\[Kappa]=a*(\[Epsilon]1*1000)^3+b*(\[Epsilon]1*1000)^2+c*\[Epsilon]1*1000+d ;
\[Sigma]sx1=If[hyp[[i,1]]=="E",Es*\[Epsilon]x,fyx1-\[Kappa]*(Acx1/Asx1)*\[Alpha]1*fctm/(1+Sqrt[500*\[Epsilon]x])];
\[Sigma]sx2=If[\[Alpha]2==0,0,If[hyp[[i,2]]=="E",Es*\[Epsilon]x,fyx2-\[Kappa]*(Acx2/Asx2)*\[Alpha]2*fctm/(1+Sqrt[500*\[Epsilon]x])]];
\[Sigma]st=If[hyp[[i,3]]=="E",Es*\[Epsilon]t,fyt-\[Kappa]*(Act/Ast)*\[Alpha]t*fctm/(1+Sqrt[500*\[Epsilon]t])];
Eq1[\[Epsilon]1_,\[Theta]_]:=Evaluate[Asx1*\[Sigma]sx1+Asx2*\[Sigma]sx2-V/Tan[\[Theta]]+\[Sigma]1*bw*z];
Eq2[\[Epsilon]1_,\[Theta]_]:=Evaluate[\[Sigma]st*Ast-(\[Sigma]2*Sin[\[Theta]]^2-\[Sigma]1*Cos[\[Theta]]^2)*bw*s];
Quiet[Sol=Check[FindRoot[{Eq1[\[Epsilon]1,\[Theta]],Eq2[\[Epsilon]1,\[Theta]]},{\[Epsilon]1,\[Epsilon]ctm+(Min[\[Epsilon]yx1,\[Epsilon]yx2]+\[Epsilon]yt-\[Epsilon]ctm)*(j-1)/(epsilon1seeds-1)},{\[Theta],seedvector[[i]]}, MaxIterations->10000],0]];
If[Length[Sol]!=0,
\[Epsilon]1=\[Epsilon]1/.Sol[[1]];
\[Kappa]=a*(\[Epsilon]1*1000)^3+b*(\[Epsilon]1*1000)^2+c*\[Epsilon]1*1000+d ;
If[\[Epsilon]1<\[Epsilon]ctm || \[Kappa]>\[Kappa]lim2, 
(* Print["            ",specimens[[spec]]," - ",hyp[[i]],"  ->  Soluci\[OAcute]n fuera de rango  \[Epsilon]1=" ,\[Epsilon]1," < \[Epsilon]ctm=", \[Epsilon]ctm,"  \[Theta]=",\[Theta]/.Sol[[2]],"  \[Kappa]=", a*(\[Epsilon]1*1000)^3+b*(\[Epsilon]1*1000)^2+c*\[Epsilon]1*1000+d ]; *)
Sol=0;
];
];
If[Length[Sol]!=0,
\[Theta]=\[Theta]/.Sol[[2]];
\[Sigma]1=\[Alpha]*fctm/(1+Sqrt[500*\[Epsilon]1]);
\[Sigma]2=(Tan[\[Theta]]+1/Tan[\[Theta]])*V/(z*bw)-\[Sigma]1;
f2max=Min[fc,fc/(0.8+170*\[Epsilon]1)];
\[Epsilon]2=\[Epsilon]c*(1-Sqrt[1-\[Sigma]2/f2max]);
\[Epsilon]t=(\[Epsilon]2*Tan[\[Theta]]^2+\[Epsilon]1)/(Tan[\[Theta]]^2+1);
\[Epsilon]x=\[Epsilon]1+\[Epsilon]2-\[Epsilon]t;
\[Sigma]st=If[hyp[[i,3]]=="E",Es*\[Epsilon]t,fyt-\[Kappa]*(Act/Ast)*\[Alpha]t*fctm/(1+Sqrt[500*\[Epsilon]t])];
solx1=FindRoot[Es*(\[Epsilon]yx1-x)== \[Kappa]*(Acx1/Asx1)*\[Alpha]1*fctm/(1+Sqrt[500*x]),{x,\[Epsilon]yx1}];
emaxx1=x/.solx1[[1]];
If[\[Alpha]2!=0,
solx2=FindRoot[Es*(\[Epsilon]yx2-x)== \[Kappa]*(Acx2/Asx2)*\[Alpha]2*fctm/(1+Sqrt[500*x]),{x,\[Epsilon]yx2}];
emaxx2=x/.solx2[[1]];
];
solt=FindRoot[Es*(\[Epsilon]yt-x)== \[Kappa]*(Act/Ast)*\[Alpha]t*fctm/(1+Sqrt[500*x]),{x,\[Epsilon]yt}];
emaxt=x/.solt[[1]];
specdiff=If[((\[Epsilon]x<=emaxx1 && hyp[[i,1]]=="E")||(\[Epsilon]x>=emaxx1 && hyp[[i,1]]=="P")) &&(\[Alpha]2==0 ||( \[Epsilon]x<=emaxx2 && hyp[[i,2]]=="E")||(\[Epsilon]x>=emaxx2 && hyp[[i,2]]=="P")) && ((\[Epsilon]t<=emaxt && hyp[[i,3]]=="E")||(\[Epsilon]t>=emaxt && hyp[[i,3]]=="P")), Append[specdiff,(\[Sigma]st-\[Sigma]stexp)^2],Append[specdiff,penal]];
(* Clear[\[Epsilon]1,\[Theta]];Print[Sol]; *)
(* Print["                     j=",j "  emaxx1=", emaxx1,If[\[Alpha]2!=0,"   emaxx2="<>ToString[emaxx2],""],"   emaxt=",emaxt]; *)
];
If[Length[Sol]==0,
specdiff=Append[specdiff,penal];
];
j++;
];
(* Print[specdiff]; *)
diff=Append[diff,Min[specdiff]];
];
(* Print["           Diferencias al cuadrado ",diff];*)
suma=suma+Min[diff];
];
Return[suma/(Length[specimens]*10^6)]; 
];
Coef[x_]:=Tan[\[Pi]*x-\[Pi]/2];
Coord[y_]:=ArcTan[y]/\[Pi]+1/2;
G2B[g_]:=Module[{i,b},
b={g[[1]]};
For[i=2,i<=Length[g],i++,
b=Append[b,If[g[[i]]==0,b[[i-1]],1-b[[i-1]]]];
];
Return[b];
];
B2Real[b_]:=Module[{i,r},
r=0;
For[i=Length[b],i>=1,i=i-1,
r=(r+b[[i]])/2;
];
Return[N[r]];
];
Real2B[r_,L_]:=IntegerDigits[Floor[r*2^L],2,L];
B2G[b_]:=Module[{i,g},
g={b[[1]]};
For[i=2,i<=Length[b],i++,
g=Append[g,If[b[[i]]==b[[i-1]],0,1]];
];
Return[g];
];


GA[n_,pmut_,pco_,repl_,G_]:=Module[{P,size,i,fit,sum},
size=80;
For[i=1,i<= G,i++,
Write[messages,"Generation "<>ToString[i]];
P= Population[n,size];
Write[messages,"   Selection ..."];
P=FPSelection[P];
Write[messages,"   Cross-over ..."];
P = Crossover[P,pco];
Write[messages,"   Mutation ..."];
P = Mutation[P,pmut];
Export[outputxls,P]
(* Export[initpop,Take[P,IntegerPart[repl*n]]]; *)
];
fit={};
For[i=1,i<=Length[P],i=i+1,
fit=Append[fit,Fitness[P[[i]]]];
];
sum=Sum[fit[[i]],{i,Length[P]}];
Write[messages,"Fitnesses: "<>ToString[fit]<>"   Mean Fitness: "<>ToString[sum/Length[P]]<>" ("<>ToString[sum*100/Length[P]]<>"%)   Max. Fitness: "<>ToString[ Max[fit]]<>" ("<>ToString[Max[fit]*100]<>"%)"];
meanfit=Append[meanfit,sum*100/Length[P]];
maxfit=Append[maxfit,Max[fit]*100];
Write[messages," "];
Write[messages,"MaxFitness: "<>ToString[maxfit]<>"    MeanFitness: "<>ToString[meanfit]];
Return[P];
];


maxfit={};
meanfit={};
P=GA[100,0.75,0.9,0.5,1];
(*P=GA[100,0.75,0.9,0.5,200];*)
Export[outputxls,P,"CSV"];
Close[messages];
